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[Document Name] Specification 

[Title of the Invention] Brain Current Source Estimating Method, Brain Current 
Source Estimating Program, and Brain Current Source Estimating Apparatus 
[Scope of Claims for Patent] 

[Claim 1] A brain current source estimating method for estimating, based on 
an electromagnetic field observed outside a scalp, a position of a current source as a 
source of said electromagnetic field existing in the brain, comprising the steps of: 

setting, in the brain, a plurality of virtual curved surfaces having depths from 
brain surface different from each other and shapes not intersecting with each other, and 
setting lattice points on each of said virtual curved surfaces; 

estimating, on each of said virtual curved surfaces, a current distribution for 
recovering said observed electromagnetic field; and 

based on an expansion of the current distribution estimated on said virtual 
curved surfaces and a difference between the electromagnetic field recovered based on 
said current distribution and said observed electromagnetic field, identifying a virtual 
curved surface at which said expansion and said difference attain relative minimums 
among said plurality of virtual curved surfaces as a true curved surface on which said 
current source exists. 

[Claim 2] The brain current source estimating method according to claim 1, 
wherein 

said step of estimating a current distribution includes the step of 

determining posterior probability by Bayesian estimation method from prior 
distribution and observation data of said electromagnetic field, and 

said step of identifying as a true curved surface on which said current source 
exists includes the step of 

identifying a virtual curved surface of which corresponding said posterior 
probability attains the maximum, among said virtual curved surfaces. 

[Claim 3] The brain current source estimating method according to claim 2, 

- 1 - 



JP200 1-4005 19 



wherein 

said step of estimating a current distribution includes the step of 

identifying a first virtual curved surface closest to said brain surface and having 
posterior probability attaining a relative maximum, among said plurality of virtual 
curved surfaces, while successively moving from a virtual curved surface on the side of 
the brain surface to a deeper side, and 

said step of identifying a curved surface as a true curved surface on which said 
current source exists includes the steps of 

identifying a localized current distribution corresponding to a point of relative 
maximum of said current distribution, on said first virtual curved surface, 

separating a plurality of local surfaces each including said localized current 
distribution, and 

fixing, among said plurality of local surfaces, local surfaces other than a local 
surface as an object of identification, moving said local surface as an object of 
identification in the depth direction, and identifying positions where said posterior 
probability attains the relative maximum, successively from the side closer to said brain 
surface. 

[Claim 4] The brain current source estimating method according to claim 3, 
wherein 

in said step of estimating a current distribution, said current distribution is 
estimated with a first spatial resolution, 

said method further comprising the step of 

re-estimating said current distribution with a second spatial resolution higher 
than said first spatial resolution and resolution of said plurality of virtual curved 
surfaces in the depth direction being improved. 

[Claim 5] A program for a computer for estimating, based on an 
electromagnetic field observed outside a scalp, a position of a current source as a source 
of said electromagnetic field existing in the brain, to have the computer execute the 
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steps of: 

setting, in the brain, a plurality of virtual curved surfaces having depths from 
brain surface different from each other and shapes not intersecting with each other, and 
setting lattice points on each of said virtual curved surfaces; 

estimating, on each of said virtual curved surfaces, a current distribution for 
recovering said observed electromagnetic field; and 

based on an expansion of the current distribution estimated on said virtual 
curved surfaces and a difference between the electromagnetic field recovered based on 
said current distribution and said observed electromagnetic field, identifying a virtual 
curved surface at which said expansion and said difference attain relative minimums 
among said plurality of virtual curved surfaces as a true curved surface on which said 
current source exists. 

[Claim 6] The program according to claim 5, wherein 
said step of estimating a current distribution includes the step of 
determining posterior probability by Bayesian estimation method from prior 
distribution and observation data of said electromagnetic field, and 

said step of identifying as a true curved surface on which said current source 
exists includes the step of 

identifying a virtual curved surface of which corresponding said posterior 
probability attains the maximum, among said virtual curved surfaces. 
[Claim 7] The program according to claim 6, wherein 
said step of estimating a current distribution includes the step of 
identifying a first virtual curved surface closest to said brain surface and having 
posterior probability attaining a relative maximum, among said plurality of virtual 
curved surfaces, while successively moving from a virtual curved surface on the side of 
the brain surface to a deeper side, and 

said step of identifying a curved surface as a true curved surface on which said 
current source exists includes the steps of 
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identifying a localized current distribution corresponding to a point of relative 
maximum of said current distribution, on said first virtual curved surface, 

separating a plurality of local surfaces each including said localized current 
distribution, and 

fixing, among said plurality of local surfaces, local surfaces other than a local 
surface as an object of identification, moving said local surface as an object of 
identification in the depth direction, and identifying positions where said posterior 
probability attains the relative maximum, successively from the side closer to said brain 
surface. 

[Claim 8] The program according to claim 7, wherein 
in said step of estimating a current distribution, said current distribution is 
estimated with a first spatial resolution, 

said program further comprising the step of 

re-estimating said current distribution with a second spatial resolution higher 
than said first spatial resolution and resolution of said plurality of virtual curved 
surfaces in the depth direction being improved. 

[Claim 9] A brain current source estimating apparatus for estimating, based on 
an electromagnetic field observed outside a scalp, a position of a current source as a 
source of said electromagnetic field existing in the brain, comprising: 

virtual curved surface setting means for setting, in the brain, a plurality of 
virtual curved surfaces having depths from brain surface different from each other and 
shapes not intersecting with each other, and setting lattice points on each of said virtual 
curved surfaces; 

current distribution estimating means for estimating, on each of said virtual 
curved surfaces, a current distribution for recovering said observed electromagnetic 
field; and 

current source identifying means for identifying, based on an expansion of the 
current distribution estimated on said virtual curved surfaces and a difference between 
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the electromagnetic field recovered based on said current distribution and said observed 
electromagnetic field, a virtual curved surface at which said expansion and said 
difference attain relative minimums among said plurality of virtual curved surfaces as a 
true curved surface on which said current source exists. 

[Claim 10] The brain current source estimating apparatus according to claim 9, 
wherein 

said current distribution estimating means includes 

posterior probability determining means for determining posterior probability by 
Bayesian estimation method from prior distribution and observation data of said 
electromagnetic field, and 

said current source identifying means includes 

virtual curved surface identifying means for identifying a virtual curved surface 
of which corresponding said posterior probability attains the maximum, among said 
virtual curved surfaces. 

[Claim 11] The brain current source estimating apparatus according to claim 
1 0, wherein 

said current distribution estimating means includes 

shallowest virtual curved surface identifying means for identifying a first virtual 
curved surface closest to said brain surface and having posterior probability attaining a 
relative maximum, among said plurality of virtual curved surfaces, while successively 
moving from a virtual curved surface on the side of the brain surface to a deeper side, 
and 

said current source identifying means includes 

localized current distribution identifying means for identifying a localized 
current distribution corresponding to a point of relative maximum of said current 
distribution, on said first virtual curved surface, 

local surface extracting means for separating a plurality of local surfaces each 
including said localized current distribution, and 
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local surface position identifying means for fixing, among said plurality of local 
surfaces, local surfaces other than a local surface serving as an object of identification, 
moving said local surface as an object of identification in the depth direction, and 
identifying positions where said posterior probability attains the relative maximum, 
successively from the side closer to said brain surface. 

[Claim 12] The brain current source estimating apparatus according to claim 
11, wherein 

said current distribution estimating means initially estimates said current 
distribution with a first spatial resolution and thereafter re-estimates said current 
distribution with a second spatial resolution higher than said first spatial resolution and 
resolution of said plurality of virtual curved surfaces in the depth direction being 
improved. 

[Detailed Description of the Invention] 

[Technical Field to Which the Invention Belongs] 

The present invention relates to a method of estimating position or distribution 
of a brain current source or sources generating magnetic fields or electric fields outside 
a scalp in correspondence with brain activities, as well as to a brain current source 
estimating program and to a brain current source estimating apparatus. 

[Prior Art] 

Thanks to remarkable development in the field of biomedical measurement 
technique, precision in measuring weak electric field (brain wave) or weak magnetic 
field (brain magnetic wave) generated from the brain, of which measurement has been 
difficult and error-prone, has been improved year by year. 

Specifically, receiving external stimuli, neural cells in the brain generate a 
current. The current results in the aforementioned weak electric field or weak 
magnetic field. Here, "brain wave" refers to the electric field generated from the brain 
by the current from the neural cells. An "electroencephalogram: EEG" refers to a 
method of measuring the brain wave. 
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The "brain magnetic wave" refers to a magnetic field generated from the brain 
by the current from the neural cells. A "magnetoencephalography: MEG" refers to a 
method of measuring the brain magnetic wave. A crucial advantage of the 
magnetoencephalography is that, as the magnetic field is almost free of any influence of 
volume conductor, it is expected that relatively accurate three-dimensional estimation 
of the position of a brain current source can be attained by measuring magnetism from 
outside one's scalp. 

In the analysis of the brain magnetic wave, an active portion of the brain is 
estimated in a non-invasive manner, by measuring the generated magnetic field from 
outside the brain. 

The magnetic field, however, is so weak that it is very susceptible to the 
influence of external magnetic field such as terrestrial magnetism. Therefore, the 
weak magnetic field is measured by a Superconducting QUantum Interference Device 
(SQUID), which is a measuring device utilizing superconductivity, within a shield that 
shuts out any external magnetic field. 

[Problems to be Solved by the Invention] 

It is noted, however, that in the field of studying algorithms for "estimating 
positions of brain current source," decisive method is non-existent at present, though 
various and many variations of initial models have been tried. 

By way of example, "dipole estimation method" as one algorithm for 
"estimating positions of brain current sources," is disclosed in Reference 1 : J. C. 
Mosher, P. S. Lewis and R. ML Leahy, IEEE Trans. Biomed. Engng. <1992> vol. 39, 
pp. 541-557. In the "dipole estimation method," however, the position of a dipole is 
estimated from observed magnetic field, assuming that the current source in the brain 
can be represented by one or a number of current dipoles, and this method is 
disadvantageous in that it is difficult to determine the number of dipoles. 

As another algorithm, "spatial filtering method" is disclosed in Reference 2: K. 
Toyama, K. Yoshikawa, Y. Yoshida, Y. Kondo, S. Tomita, Y. Takanashi, Y. Ejima, 
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and S. Yoshizawa, Neuroscience, 1999, 91 (2), pp. 405-415. In the "spatial filtering 
method," location of a brain current source is restricted in consideration of 
physiological findings, and distribution of dipoles is estimated. This method is 
disadvantageous in that accurate estimation of the depth of the current source is not 
possible. 

The present invention was made to solve the above-described problems, and an 
object of the present invention is to provide a brain current source estimating method 
that enables estimation of a position, depth direction inclusive, of a brain current source 
from observed data. 

Another object of the present invention is to provide a brain current source 
estimating method that enables estimation of positions of a plurality of brain current 
sources from observed data, even when there are a plurality of brain current sources. 

A still further object of the present invention is to provide a brain current source 
estimating program that enables estimation of a position, depth direction inclusive, of a 
brain current source from observed data. 

A still further object of the present invention is to provide a brain current source 
estimating program that enables estimation of positions of a plurality of brain current 
sources from observed data, even when there are a plurality of brain current sources. 

A still further object of the present invention is to provide a brain current source 
estimating apparatus that enables estimation of a position, depth direction inclusive, of 
a brain current source from observed data. 

A still further object of the present invention is to provide a brain current source 
estimating apparatus that enables estimation of positions of a plurality of brain current 
sources from observed data, even when there are a plurality of brain current sources. 

[Means for Solving the Problems] 

In order to attain such objects, a brain current source estimating method 
according to the present invention is a method for estimating, based on an 
electromagnetic field observed outside a scalp, a position of a current source as a source 
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of the electromagnetic field existing in the brain, including the steps of: setting, in the 
brain, a plurality of virtual curved surfaces having depths from brain surface different 
from each other and shapes not intersecting with each other and setting lattice points on 
each of the virtual curved surfaces; estimating, on each of the virtual curved surfaces, a 
current distribution for recovering the observed electromagnetic field; and based on an 
expansion of the current distribution estimated on the virtual curved surfaces and a 
difference between the electromagnetic field recovered based on the current distribution 
and the observed electromagnetic field, identifying a virtual curved surface at which the 
expansion and the difference attain relative minimums among the plurality of virtual 
curved surfaces as a true curved surface on which the current source exists. 

Preferably, the step of estimating a current distribution includes the step of 
determining posterior probability by Bayesian estimation method from prior 
distribution and observation data of the electromagnetic field, and the step of 
identifying as a true curved surface on which the current source exists includes the step 
of identifying a virtual curved surface of which corresponding posterior probability 
attains the highest, among the virtual curved surfaces. 

Preferably, the step of estimating a current distribution includes the step of 
identifying a first virtual curved surface closest to the brain surface and having posterior 
probability attaining a relative maximum, among the plurality of virtual curved surfaces, 
while successively moving from a virtual curved surface on the side of the brain surface 
to a deeper side, and the step of identifying a curved surface as a true curved surface on 
which the current source exists includes the steps of identifying a localized current 
distribution corresponding to a point of relative maximum of the current distribution, 
on the first virtual curved surface, separating a plurality of local surfaces each including 
the localized current distribution, and fixing, among the plurality of local surfaces, local 
surfaces other than a local surface as an object of identification, moving the local 
surface as an object of identification in the depth direction, and identifying positions 
where the posterior probability attains the relative maximum, successively from the side 
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closer to the brain surface. 

Preferably, in the step of estimating a current distribution, the current 
distribution is estimated with a first spatial resolution, and the method further includes 
the step of re-estimating the current distribution with a second spatial resolution higher 
than the first spatial resolution and resolution of the plurality of virtual curved surfaces 
in the depth direction being improved. 

According to another aspect of the present invention, a program for a computer 
for estimating, based on an electromagnetic field observed outside a scalp, a position of 
a current source as a source of the electromagnetic field existing in the brain, has the 
computer execute the steps of: setting, in the brain, a plurality of virtual curved surfaces 
having depths from brain surface different from each other and shapes not intersecting 
with each other and setting lattice points on each of the virtual curved surfaces; 
estimating, on each of the virtual curved surfaces, a current distribution for recovering 
the observed electromagnetic field; and based on an expansion of the current 
distribution estimated on the virtual curved surfaces and a difference between the 
electromagnetic field recovered based on the current distribution and the observed 
electromagnetic field, identifying a virtual curved surface at which the expansion and 
the difference attain relative minimums among the plurality of virtual curved surfaces 
as a true curved surface on which the current source exists. 

Preferably, the step of estimating a current distribution includes the step of 
determining posterior probability by Bayesian estimation method from prior 
distribution and observation data of the electromagnetic field, and the step of 
identifying as a true curved surface on which the current source exists includes the step 
of identifying a virtual curved surface of which corresponding posterior probability 
attains the highest, among the virtual curved surfaces. 

Preferably, the step of estimating a current distribution includes the step of 
identifying a first virtual curved surface closest to the brain surface and having posterior 
probability attaining a relative maximum, among the plurality of virtual curved surfaces, 
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while successively moving from a virtual curved surface on the side of the brain surface 
to a deeper side, and the step of identifying a curved surface as a true curved surface on 
which the current source exists includes the steps of identifying a localized current 
distribution corresponding to a point of relative maximum of the current distribution, 
on the first virtual curved surface, separating a plurality of local surfaces each including 
the localized current distribution, and fixing, among the plurality of local surfaces, local 
surfaces other than a local surface as an object of identification, moving the local 
surface as an object of identification in the depth direction, and identifying positions 
where the posterior probability attains the relative maximum, successively from the side 
closer to the brain surface. 

Preferably, in the step of estimating a current distribution, the current 
distribution is estimated with a first spatial resolution, and the program further includes 
the step of re-estimating the current distribution with a second spatial resolution higher 
than the first spatial resolution and resolution of the plurality of virtual curved surfaces 
in the depth direction being improved. 

According to a still further aspect of the present invention, a brain current 
source estimating apparatus for estimating, based on an electromagnetic field observed 
outside a scalp, a position of a current source as a source of the electromagnetic field 
existing in the brain, includes: virtual curved surface setting means for setting, in the 
brain, a plurality of virtual curved surfaces having depths from brain surface different 
from each other and shapes not intersecting with each other and setting lattice points on 
each of the virtual curved surfaces; current distribution estimating means for estimating, 
on each of the virtual curved surfaces, a current distribution for recovering the observed 
electromagnetic field; and current source identifying means for identifying, based on an 
expansion of the current distribution estimated on the virtual curved surfaces and a 
difference between the electromagnetic field recovered based on the current distribution 
and the observed electromagnetic field, a virtual curved surface at which the expansion 
and the difference attain relative minimums among the plurality of virtual curved 
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surfaces as a true curved surface on which the current source exists. 

Preferably, the current distribution estimating means includes posterior 
probability determining means for determining posterior probability by Bayesian 
estimation method from prior distribution and observation data of the electromagnetic 
field, and the current source identifying means includes virtual curved surface 
identifying means for identifying a virtual curved surface of which corresponding 
posterior probability attains the highest, among the virtual curved surfaces. 

Preferably, the current distribution estimating means includes shallowest virtual 
curved surface identifying means for identifying a first virtual curved surface closest to 
the brain surface and having posterior probability attaining a relative maximum, among 
the plurality of virtual curved surfaces, while successively moving from a virtual curved 
surface on the side of the brain surface to a deeper side, and the current source 
identifying means includes localized current distribution identifying means for 
identifying a localized current distribution corresponding to a point of relative 
maximum of the current distribution, on the first virtual curved surface, local surface 
extracting means for separating a plurality of local surfaces each including the localized 
current distribution, and local surface position identifying means for fixing, among the 
plurality of local surfaces, local surfaces other than a local surface as an object of 
identification, moving the local surface as an object of identification in the depth 
direction, and identifying positions where the posterior probability attains the relative 
maximum, successively from the side closer to the brain surface. 

Preferably, the current distribution estimating means initially estimates the 
current distribution with a first spatial resolution and thereafter re-estimates the current 
distribution with a second spatial resolution higher than the first spatial resolution and 
resolution of the plurality of virtual curved surfaces in the depth direction being 
improved. 
[Embodiments] 

Embodiments of the present invention will be described hereinafter with 
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reference to the figures. 

As already described, magnetoencephaography: MEG and 
electroencephalogram: EEG have been known as methods of monitoring brain activities 
from the outside. In the following, description will be given mainly focusing on MEG. 
It is noted, however, that the present invention is also applicable to EEG, when the 
magnetic field in MEG is replaced by the electric field. 

The term "electromagnetic field" originally refers to co-existence of "electric 
field" and "magnetic field." In the present specification, however, an expression 
"observe an electromagnetic field" will be used generally, where the physical amount 
"to be observed" is an "electric field" or a "magnetic field." 

Fig. 1 is a schematic illustration of an exemplary configuration of an MEG 

system. 

MEG 12 includes an array of multi-channel SQUID flux meter (super sensitive 
magnetometer), for measuring a magnetic field on the scalp surface of a subject 10. 
Receiving the result of measurement by MEG 12, a computer system 20 analyzes the 
result of measurement, and performs a process for estimating the position of a brain 
current source. 

The present invention relates to software processing by computer system 20. It 
is noted, however, that part of or all of the process may be performed by hardware to 
increase the speed of processing. 

[Principle of brain current source estimation] 

Prior to detailed description of the "brain current source estimating method" in 
accordance with the present invention, the principle and outline of the estimation 
method of the present invention will be summarized. 

(Recovery of electromagnetic field by current distribution on a curved surface) 
It is well-known as the principle of electromagnetic shield that when a current 
source is completely surrounded by a metal surface, electromagnetic field ceases to 
exist outside the metal surface. When viewed from a different point, this means that 
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an electromagnetic field identical to that formed by the current source can be formed 
outside the metal surface by causing an appropriate current flow over the metal surface. 
On the contrary, the electromagnetic field formed by the current source located outside 
the metal surface cannot be fully recovered, whatever current is caused to flow over the 
metal surface. This fact will be referred to as the "principle of electromagnetic field 
recovery." 

Fig. 2 is a schematic illustration representing an electromagnetic field generated 
by a current source, observed on an appropriate curved surface. 

As can be seen from Fig. 2, when a virtual curved surface is prepared between 
an observing surface and the current source, it is possible to recover the electromagnetic 
field formed by the current source on the observing surface by causing an appropriate 
flow of current on the virtual curved surface, in accordance with the "principle of 
electromagnetic field recovery" described above. 

Fig. 3 is a schematic illustration representing a relation between a current source 
in the brain and a plurality of virtual curved surfaces. 

Referring to Fig. 3, considering that the electromagnetic field formed by the 
current attenuates in inverse proportion to the square of distance, the expansion of 
current on the virtual curved surface (virtual curved surface 2) equivalent to the current 
source becomes wider as the virtual curved surface is further away from the current 
source. Therefore, the expansion of the current on virtual curved surface 1 is wider 
than that on virtual curved surface 2. 

On a virtual curved surface (virtual curved surface 3) on the side opposite to the 
observing curved surface with respect to the current source, it is impossible to fully 
recover the electromagnetic field formed on the observing surface by the current source. 

According to the present invention, based on the principle described above, the 
current source is estimated from the observed data of the electromagnetic field on the 
observing surface. 

(Principle of current source estimation) 
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Assume that a magnetic field (or an electric field) formed by a current source 
generated in the brain is observed on an observing surface that is close to the surface of 
the brain, as shown in Fig. 2. 

A virtual curved surface in the brain is considered, and current distribution on 
the virtual curved surface that recovers the observed magnetic field is calculated. 
When the virtual curved surface is moved from the surface to the inner side of the brain 
with the radius made gradually smaller, the expansion of the current distribution 
recovering the observed magnetic field becomes smaller, and it becomes the smallest 
when the virtual surface encompasses the true current source. When the virtual curved 
surface is further moved to be deeper than the current source, the expansion of current 
distribution on the virtual curved surface comes to be wider again, and the difference 
between the magnetic field generated by the current and the observed magnetic field 
also becomes larger. 

Accordingly, the depth of the current source can be identified by reviewing the 
expansion of the current distribution recovering the observed magnetic field and the 
error in recovery of the magnetic field. Further, by calculating the current distribution 
on the virtual curved surface at the depth identified in this manner, the expansion of the 
current source can also be found. The foregoing is the principle of the present 
invention. 

(Identification of current source depth by Bayesian estimation) 
In order to specifically implement the principle of current source estimation 
described above, in the present invention, a procedure based on the following 
probability theory is employed. The procedure will be summarized in the following. 

What can be observed by MEG or EEG is several tens of magnetic fields 
(electric fields) existing near the surface of the brain. In order to approximate the 
current distribution on the virtual curved surface, lattice points are set on the virtual 
curved surface, and a current dipole (or an appropriate current source model) is 
allocated to each lattice point. In order to estimate the current distribution with high 
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resolution, it is necessary to increase the number of lattice points to increase the density 
of the lattice points. 

When the number of lattice points on the virtual curved surface is increased to 
be larger than the number of observation points, a unique solution cannot be determined. 
When estimation points is larger in number than the observation points, the number of 
parameters to be estimated becomes larger than the number of observed data, and hence, 
the observed magnetic field would be better recovered even on a virtual curved surface 
that is positioned deeper than the current source. 

In order to cope with such problems, in the present invention, the current 
distribution on the virtual curved surface is estimated utilizing Bayesian estimation 
theory. At the time of Bayesian estimation, prior distribution that represents prior 
information of the current source is used. Specifically, a prior distribution that leads 
to the smallest possible expansion of current distribution is introduced. Namely, a 
prior distribution is introduced in which a current dipole on a lattice point of which 
magnitude cannot be well identified only from the observed data comes to have a 
magnitude close to zero. This can be realized by introducing hierarchical prior 
distribution referred to as "Automatic Relevance Determination: ARD" prior 
distribution. 

It is impossible, however, to analytically calculate posterior probability 
distribution from the ARD prior distribution and observed data. Therefore, in the 
present invention, variational Bayes method is used as a method of calculating posterior 
probability distribution by approximation, as will be described later. It is noted that 
other method of approximation such as Monte Carlo method may be used. 

By Bayesian estimation using ARD prior distribution, it becomes possible to 
obtain a current distribution on the virtual curved surface that recovers the observed 
data and has smallest possible expansion. Further, by comparing model posterior 
probabilities when estimations are made using virtual curved surfaces of different 
depths, the depth of the current source can be identified. 
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The logarithm of the model posterior probability can be represented as a sum of 
log-likelihood term and model complexity term. The log-likelihood term becomes 
larger as recovery error becomes smaller. 

The model complexity term becomes larger when the number of effective 
current dipoles (that is, having a magnitude not smaller than an appropriate threshold) 
on the lattice points becomes smaller. As already described, the recovery error and the 
expansion of the current distribution on the virtual curved surface (number of effective 
current dipoles) become the smallest when the depth of the virtual curved surface 
matches the current source. Thus, it can be understood that the model posterior 
probability becomes the largest at this time. In other words, it can be understood that 
the current source exists on the virtual curved surface at a depth at which the model 
posterior probability becomes the largest. 

In summary, the present invention provides a method of estimating the position 
of the brain current source, depth direction inclusive, from the observed data of MEG 
and EEG. 

The basic idea of the present invention comes from the fact that an 
electromagnetic field generated by a current source can be recovered by causing an 
appropriate current flow over a curved surface existing between the current source and 
an observing surface, and that the expansion of current distribution on the curved 
surface becomes smaller as the curved surface comes closer to the current source, as 
described above. The present invention is characterized in that, based on this idea, the 
position of the current source including the depth direction is estimated by considering 
the fact that the model posterior probability becomes the largest when the curved 
surface encompasses the current source in Bayesian estimation of the current 
distribution on the curved surface that recovers the observed data, that is, by 
considering the model posterior probability. 

(Where there are a plurality of current sources) 

Though description has been made on one current source, the method is also 
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applicable even when there are a plurality of current sources. 

The electromagnetic field generated by a current attenuates in inverse proportion 
to the square of distance, and therefore, the current source closest to the brain surface 
has the largest influence on the observed magnetic field on the brain surface. 
Therefore, it is possible to identify the current sources one by one in order, starting 
from the one closest to the brain surface. 

When the virtual curved surface is moved deeper from the brain surface, the 
model posterior probability attains the relative maximum near a current source closest 
to the brain surface (which will be referred to as the first current source). When there 
are two or more current sources, there will be a plurality of local sets of current dipoles 
corresponding in number to the current sources, in the current distribution on the virtual 
curved surface. 

The local set of current dipoles is referred to as localized current distribution. 
A local surface including individual localized current distribution is separated and 
moved in the depth direction to find the model posterior probability. When the local 
surface that corresponds to the first current source is moved, the model posterior 
probability attains to the relative maximum at the depth of the first current source. 
When other local surfaces are moved, however, the model posterior probability never 
attains to the relative maximum at the depth of the first current. Thus, the position of 
the first current source, that is, the position of the current source closest to the brain 
surface can be identified. 

In order to find the second deepest current source, the local surface 
corresponding to the first current source is fixed, and the depths of the remaining local 
surfaces are aligned and gradually made deeper. Then the model posterior probability 
attains the relative maximum at the second deepest current source. When the 
individual local surface is moved in the depth direction again, the model posterior 
probability attains to the relative maximum at the depth of the second current source, 
only when the local surface corresponding to the second current source is moved. In 
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this manner, the position of the second current source can be identified. The third and 
the following current sources can also be identified in the similar manner. 

The method is advantageous over the method in which the depth of individual 
local surface is moved independently, in that the time for computation can significantly 
be reduced. 

(Method in which resolution is increased gradually) 

According to the method of estimating brain current source of the present 
invention, it is possible to identify position of each of the current sources starting from 
the one closest to the brain surface in the above described manner, and in addition, it is 
possible to gradually increase the resolution of current source estimation. 

First, a position of a current source is roughly estimated with a low resolution 
(with a small number of lattice points on the virtual curved surface and a small number 
of sample points in the depth direction). In this stage, the position of a local surface 
corresponding to each current source is approximately determined. 

Next, estimation is made with higher resolution. At this time, the area of the 
local surface has become smaller as compared with the original virtual curved surface, 
and hence the resolution is naturally higher when the same number of lattice points is 
used. Accuracy in estimating the current source position can be improved by 
increasing the resolution in the depth direction as well. If the current distribution is 
further localized when the resolution is made higher, it is possible to estimate again 
with local surface made smaller and the resolution made still higher. 

On the contrary, if the expansion of current distribution does not much vary 
even when the resolution is improved, it means that the current source expands wide. 
In this manner, the resolution can be adjusted in accordance with how the current 
source expands. 

[Specific procedure of current source estimation] 

In the following, specific procedure for identifying the position of a "brain 
current source" will be described in detail, in accordance with the outline above. 
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[(I) Preparation for current source estimation] 

First, preparation for the process procedure for estimating the current source 
will be described. 

Specifically, for the estimation of current source, the following procedures must 
be taken in advance. 

(1-1) Determination of the shape of virtual curved surface and current model 

(1-2) Further, in order to estimate the virtual current distribution while moving 
the virtual curved surface, it is necessary to determine sample points on the virtual 
curved surface and sample points in the depth direction. 

(1-3) Further, as will be described in detail later, it is necessary to determine in 
advance parameters to designate the distribution shape of hierarchical prior distribution 
for estimating current distribution as initial values. 

In the following, the aforementioned procedure of (1-1) Determination of the 
shape of virtual curved surface and current model will be described in grater detail. 

(1-1-1) Determination of the shape of virtual curved surface 

The simplest shape of the virtual curved surface is obtained by regarding the 
brain as a hemisphere and assuming various hemispheres of different radii to be the 
virtual curved surfaces. 

When a location where existence of a brain current is highly likely such as the 
cerebral cortex has been known in advance by other measuring methods such as 
Magnetic Resonance Imaging: MRI, the shape of the virtual curved surface may be 
determined based on such information. 

In this case, the shapes of virtual curved surfaces having different depths may 
generally differ. It is necessary, however, to determine the shapes not to intersect with 
each other. 

(1-1-2) Determination of current model 

As a current model on the virtual curved surface, let us consider a current dipole 
model. It is also possible to consider other current models. 
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Consider appropriate lattice points (sample points) {X n |n=l, N} on the 

virtual curved surface. Assume a current dipole of which current intensity is j n on each 

lattice point. Here, the magnetic field formed by the current dipole j n at a lattice point 

X n on an observation point Yj (i=l, I) on the brain surface is given by the following 

expression. 

Lx^-XJ 
|Y,-X n | 3 

Here, (a represents magnetic permeability, and by way of example, for a vector 
X n , the expression |X n | represents the absolute value of vector X n . 

Here, the magnetic field formed by the current dipole {j n |n=l, N} of the 
virtual curved surface on an observation point Yi (i=l, I) is represented by the 
following equation. 

|Y,- X „r 

Assuming that the direction of the magnetic field observed at the observation 
point Yj is a vector Si, c (c=l, C), component Bi, c in the direction of vector Si, 0 of the 
magnetic field at this point can be given as 



B *" 2> |Y,-X.f 



n=l 



Further, when a local gradient of a magnetic field is to be measured as a 
magnetic field to be observed, a differentiation of the equation above by Yj will be 
observed. 

When the direction of a current dipole at a lattice point X n (position vector) is 
restricted, the current dipole can be represented by the following equation, with a 
possible independent direction of the current dipole being b n ,d (d=l, D). In this 
equation, a case where D = 3 corresponds to a case where there is no restriction on the 
direction. 
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D 

Jn = 5j jn,d * kn,d 

d=l 

In summary, a magnetic field formed by the current dipole at a lattice point {X n 
|n=l, N} of the virtual curved surface on the observation point Yj can be given as 

Bi, c =Z Zjn, rf G(i,c;n,d) 
G(i,c;n,d) = n j— — — ja ~ 

| Y i ~ X n| 

Here, j n ,d is an independent component of the current dipole at the lattice point 
X n . Further, the function G (i, c; n, d) represents the component in Si, c direction of the 
magnetic field formed by the current dipole b n ,a at the lattice point X n on the 
observation point Yi. 

The problem of estimating current distribution may be considered as a problem 
of estimating a current distribution on the virtual curved surface {j n ,d I n=l, N; 
d=l, D} from the observed magnetic field {Bi, 0 |i=l, I; c=l, C}. 

The measured value of the magnetic field at the observation point described 
above may be given by the following matrix expression. 

B = G J 

B = (B i c |z = l, c = l,--sC): (J xC) dimensional vector 

J = (j^ d \ n = l,-~,N;d = !,--, D) :(NxD) dimensional vector 

G = (G(i,c;n,d)| i = c = 1,-.,C; n = l.-.N; d = 1,-,D) 

: (I x C) x (N x D) dimensional matrix 

(Current source probability model) 

With the current model determined in this manner, the following "current 
source probability model" is introduced for such current distribution estimation. 
(1-1-3) Current source probability model 

A probability model for the current model on the virtual curved surface 
described above will be considered as follows. It is assumed that the observed 
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magnetic field is represented as a sum of the magnetic field formed by the current 
distribution J on the virtual curved surface and the observation noise. Further, it is 
assumed that the observation noise is represented as a sum of a bias component 
irrespective of a measurement site and Gaussian noise having an independent variance 
a 2 at each measuring point. 

Specifically, the observed magnetic field is considered as 

B = G- J + § 

J — (Jvirtual curved surface? Jo) 

J virtual curvedsurface= (in,d \ ^ -\ 9 "%N; d =l y "' y D) 

Jo: Bias component of noise 

G = (G(i,c;n,d) | i = c = 1,-,C; n = 1.....N; d = 1,-,D) 

^ = (^ c |i = l s ..-,I;c = l,.-,C) 

: Gaussian noise with each component having independent variance a 1 

The probability distribution for the current model may be considered as follows. 

First, a virtual curved surface at a specific depth, or a set of local surfaces with 
the depth of each local surface identified, will be represented by M. When a current 
distribution J on the virtual curved surface M is given, the probability P (B | J, p, M) 
that the observed magnetic field is B is represented by the following equation, where p= 
l/o 2 

P( B I J, p, M) = exp [- - p|B - G • J| 2 + - ( I ■ C) log ( p / 2ti )] 

P = l/a 2 

(1-1-4) Hierarchical prior distribution 

As already described, ARD hierarchical prior distribution will be used as the 
prior distribution for the current distribution J on the virtual curved surface M. 
The prior distribution for the current distribution J before observation 
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(probability of J being realized) is represented by the following equation, under the 
assumption of ARD hierarchical prior distribution. 

P 0 (J | a, 13, M) = 

^ n-1 d=1 n=1 

or oi 0 2 +^-log(jS or 0 /27r)] 

P o (0 | r, M)=r(£ I 1 / zr , r^o) 

Here, T(...) represents gamma distribution, which is defined below. In the 
expression below, r(y 0 ) is a gamma function, which is also defined below. 

IXPIb, Yo ) - ^(YoP/b) Yo 7^e"^ /b 

n^^dt^-V' 

Further, a and t are hyper parameters introduced to model the current 
distribution J and prior distribution for inverse variance of noise p. Hierarchical prior 
distributions P 0 (a | M) and P 0 (t|M) for a and t are 

T> 0 (a\M) = f[r(a n \a 9 r a0 ) 
P 0 (x|M)-r(T|T 0> y x0 ) 

In the equations above, a value for each parameter for a and t will be described 
in detail later. 

(1-1-5) Bayesian estimation 

When data B of a magnetic field is observed, the posterior probability 
distribution P (J|B, M) that the current distribution is J can be calculated in the 
following manner, using Bayesian theorem. 

P(J|B,M) = Jdp da dx P(J, p,a,T |B,M) 

Here, the posterior probability distribution P (J, p, a, t|B, M) is given as 
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P(J, p,a,x|B,M) = 



P(B j J, P,M)P 0 ( J |a, p,M)P Q (P |t,M)P 0 (a |M)P 0 (t |M) 

P(B|M) 



P(B|M) = 

|dJdpdadTP(B | J, p,M)P 0 (J |a,p,M)P 0 (P |t,M)P 0 (a |M)P 0 (t |M) 



Using this posterior probability distribution, an expected value E[J|B, M] of 
current distribution is given as 



Further, P (B|M) represents marginal likelihood of the virtual curved surface 
model M. When the depth of a current source is estimated, a number of models 
having different depths are compared. At this time, it is assumed that there is no prior 
information about the depth. Specifically, in the following, it is assumed that P(M) = 
constant. 

When the observation data B is given, the probability that the model M is a true 
model, that is, the model posterior probability P(M|B), is in proportion to the model 
marginal likelihood P(B|M) under the assumption described above, and hence, the 
following relation holds. 

P(M|B)ocP(B|M) 

(1-1-6) Variational Bayes method 

It is generally impossible to analytically calculate the model marginal likelihood 
when the probability model and the hierarchical prior distribution are given. 

Therefore, as a method of calculating by approximation the model marginal 
likelihood, variational Bayes method is used. It is possible to use other method of 
approximation, such as Monte Carlo method and Laplacian approximation. 

The "variational Bayes method" will be briefly outlined, and specific procedure 
will be described later. 




-25 - 



JP2001-400519 



In order to calculate a true posterior distribution P (J, p, a, t|B, M) by 
approximation, a trial posterior distribution Q(J, P, a, t) is considered. 

Closeness between the two probability distributions P (J, p, a, t|B, M) and Q(J, 
P, a, t) can be calculated by using K-L distance given by the following expressions. 

KL(P||Q) 

= JdJdpdad,Q(J,p,a,,)log [?? ™^] 

= logP(B|M)-F(Q)>0 

F(Q) - JdJdpdadxQ( J, p,a,x)x 

■ r P(BlJ ? P ? M)P 0 (J|a,p > M)P Q (p|T 3 M)P 0 (a|M)P 0 (TlM) 
gL Q(J,P,a,x) 

The K-L distance attains to zero only when the two distributions are equal to 
each other, and otherwise it always assumes a positive value. 

When a free energy F(Q) for the trial posterior distribution Q is defined in the 
expression above, the following inequality results. 

F(Q)<logP(B|M) 

Specifically, the distribution Q(J, p, a, x) that maximizes the free energy F(Q) 
becomes equal to the true posterior distribution P (J, p, a, t|B, M), and at this time, the 
free energy is equal to the marginal log-likelihood. 

According to variational Bayes method, the form of function Q(J, p, a, x) is 
restricted to the following form, to maximize the free energy. 

Q(J,P,a,x)=:Q J (J,P)Q a (a,x) 

By alternately repeating the step of fixing the second term Qcc on the right side 
of the equation above and maximizing F(Q) with respect to the first term Qj on the right 
side and the step of fixing the first term Qj and maximizing F(Q) with respect to the 
second term Qcc, a distribution Q* is obtained that attains the relative maximum of free 
energy F(Q). 

[(H) Procedure of current source estimation] 
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A specific procedure for estimating the current source after the preparation 
above will be described hereinafter with reference to the figures. 

Fig. 4 is a flow chart representing an overall procedure of the brain current 
source position estimating method in accordance with the present invention. 

Referring to Fig. 4, first, when the process of estimating a position of a brain 
current source starts (step SI 00), values of parameters for designating the shape of 
hierarchical prior distribution for estimating the current distribution and an initial value 
of a variable to be estimated are determined (step SI 02). 

Thereafter, initial estimation of the current source is performed, using an initial 
resolution, so as to extract candidates of the current source (step SI 04). 

Then, among the current source candidates extracted in this manner, a position 
of a current source that is closest to the brain surface is identified (step SI 06). 
Thereafter, depths of other current sources are identified successively (step SI 08). 

After the positions of current sources are identified with the initial resolution in 
the above-described manner, re-estimation of the positions of the current sources is 
performed with the spatial resolution increased (step SI 10), and the process ends (step 
SI 12). 

Processes of respective steps of Fig. 4 will be described in greater detail in the 
following. 

(II-l) Initial estimation of current source using initial resolution (extraction of 
current source candidates) 

Fig. 5 is a flow chart representing the process of step SI 04 among the steps 
shown in Fig. 4, that is, the process of initial estimation of the current source using an 
initial resolution. 

Referring to Fig. 5, first, sample points in the depth direction with the initial 
resolution are determined to be {R k |k=l, K}. Current distribution is estimated for 
the virtual curved surface at each depth Rk. Further, based on the current distribution 
estimated in this manner, posterior probability for each depth Rk is calculated. The 
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posterior probability corresponds to the free energy value, which will be described later 
(step S202). 

For the current distribution at the depth R M at which the posterior probability 
calculated in this manner becomes the highest, a relative maximum point of current 
intensity is found (step S204). 

Further, a local surface that encompasses each relative maximum point is 
determined. Assuming that there are L local surfaces, each local surface will be the 
candidate of localized current source (step S206). 

In the following, the process of step SI 02 of Fig. 4 and the process of step S202 
of Fig. 5 that follows, will be described in greater detail. 

(EH-1) Determination of parameter values for designating shape of hierarchical 
prior distribution for estimating current distribution and initial values of variables to be 
estimated 

As described above, estimation of a current source is performed using 
variational Bayes method. Therefore, the process of step SI 02 shown in Fig. 4 is 
performed in the following manner, with the application of variational Bayes method. 

The parameters designating the shape of hierarchical prior distribution are 
determined as follows. 



(more generally, 0 < y ^ 0 < 1) 



7r0 



(more generally, 0 < y r0 < 1) 



= 0.1 



(more generally, 0 < y a0 < 1) 



T 0 =K t V ar (B), k x =1 



(more generally, K t ~l) 




a = K a • -^Tr (G' • G), K a = 0, 1 (more generally, 0 < K a < 1) 



N 
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N = D- N + l 

Further, based on the equations above, each variable is initialized in the 
following manner. 

Yx = Y t0 + Y po 

Y« = Y«o + -jD 

a n =a 
f = t 0 
S G =G G 

(II-1-2) Specific process of variational Bayes method for estimating current 
distribution 

(1) Calculation of expected values of parameters J, p (M-step process) 
Here, expected values of parameters J and P are calculated. By this process, 
the free energy F(Q) for Qj is maximized. 

By defining a diagonal matrix A as follows, Qj is derived in accordance with the 
following equations. In the following equations, an expected value of a variable is 
represented by the variable name with " — " attached thereabove. 

A(n,d;n',cT) = 5 nn ,d M ,a n («, n' = \,-,N; d,d' = l,-,D) 

A(n,d;0) = A(0;n',d') = 0 
A(0;0) = a 0 
2 = 2 G +A 

j=r'G'B 

P = Yp [i| B - G ■ J | 2 + 1 J' A J + Ypo f r 1 

Q J (J,P) = Q J (J|P)Q P (P) 
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Q J (J|P) = exp[-ip(J-J)'2(J-J) + ilog|s| + |Nlog(p/27 t )] 

Qp(P) = r(p|p, Yp ) 

(2) Calculation of expected values of hyper parameters, that is, calculation of 
expected values of parameters a, t (process for H-step) 

Following the M step, expected values of hyper parameters a and t are 
calculated. In this step, a process is performed to maximize the free energy F(Q) for 

Qc, 

The procedure can be represented as 

1 _d _ d 

<x n =ra[r aQ s- } + — Bl i n .d 2 + Kl^Xn. d:n. d)]" 1 

2 d=1 <1=1 

(n=1, N) 

2 2 

_ _ 1 N _ 

Q a (a. r)=r(r | r, r r >r(tf 0 l r ff0 + — ) II TCaJ <*„. r a ) 

2 n=1 

(3) Calculation of free energy 

Using Qj and Q a calculated through the M step and H step as described above, 
the free energy is calculated in the following manner. 
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F=LP + H J + H j8 +H a + H r 

LP = —J" I B-G-J] z — Trl" 1 I G + — (I-C) C < log $ > - log2 7T ) 

2 2 2 

<io g /3> = iogT+tf (r^)-iogr 5 

^ ( r ) = - j ELL2122-, ^ : digamma function 

dr 

1 — 1 — 

Hj== [TrZ-'A — log I Z H A| -N] $ J'AJ 

2 2 

H, = *(r*. r^o) + r j80 C<iogiS > + <io g r >- r7+i] 
<iogr >=iogT+^(r r )-iogr T 

cd < r , r 0 ) ^ Oog r ( r > - r ^ ( r ) + r ) ~ Oog r ( r o) — r O iog r 0 + r 0 ) 
H r =cp(r r , r r0 ) + r r0 C<'ogr >-iogT 0 -Tr^- 1 + i] 

N 

H ff =Z [<l>(r tf , r « 0 )+ r a0 «*°&cx n > — loga— a" n a~ 1 + i )] 

n=1 1 _ 

+ r a0 + — . r« 0 ) + r a0 (<Joga 0 >— loga— Ofoa-'-hl ) 

<logc* n >==logl7 n -{-# ( r a >— Jogr « <n==1 fi '", N) 

<iogc* 0 >-ioga" 0 +^ ( r „ 0 + — >—iog( r a0 +- — ) 

In this manner, the process from the M step to calculation of free energy 
described above is repeated until the value of free energy F converges. The value of 
free energy F(Q) after convergence gives an approximation of marginal log-likelihood 
logP(B|M). 

Further, model log-posterior probability differs only by the constant from log P 
(B|M), and therefore, the model having the maximum model posterior probability is the 
same as the model of which free energy is the highest, in accordance with the above- 
described approximation. 

By the above-described procedure, posterior probability for the depth Rk can be 
calculated. By performing the processes of steps S204 and S206 of Fig. 5 described 
above accordingly, it is possible to find candidates of current sources localized to 
respective local surfaces. 
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(II-2) Identification of current source closest to the brain surface 

Next, the process of step SI 06 of Fig. 4, that is, the process for identifying the 
current source closest to the brain surface among the candidates of current sources 
found in the manner as described above, will be described. 

Fig. 6 is a flow chart representing a process for identifying a current source 
closest to the brain surface. 

First, as the initial value, the value of a variable 1 is set to 1 (step S302). Then, 
the variable 1 is compared with a possible maximum value L of variable 1 (step S304), 
and when the variable 1 is not larger than the maximum value L, the process proceeds to 
step S3 06, and when the variable 1 exceeds the maximum value L, the process proceeds 
to step S324 (step S304). Specifically, the process from step S306 to step S322 is 
repeated from 1=1 to l^L. 

In step S3 06, depth of a local surface other than the 1th local surface is fixed at 
the depth Rm calculated in step S204 of Fig. 5. 

The value of variable k is set to 1 (step S3 08). Thereafter, the variable k is 
compared with a possible maximum value K of variable k (step S3 10), and when the 
variable k is not larger than the maximum value K, the process proceeds to step S3 12, 
and when the variable k exceeds the maximum value K, the process proceeds to step 
S320. 

In step S3 12, the depth of the 1th local surface is first determined to be Rk. 
Thereafter, current distribution is estimated for the set of L local surfaces (step 

S3 14). 

Thereafter, the posterior probability (free energy) for this arrangement of the 
local surfaces is calculated (step S3 16). By incrementing the value of variable k by 
only 1 (step S3 18), the process returns to step S3 10. 

The process is performed for each of the local surfaces having the depths from 
k=l to k=K, and then, the depth Rm(1) of the 1th local surface that provides the highest 
posterior probability is calculated (step S3 20). By incrementing the value of variable 1 
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only by 1 (step S3 22), the process returns to step S3 04. 

In this manner, among the depths RmG) of the local surfaces calculated for each 
variable 1, 1 closest to the brain surface (shallow) is detected, which is denoted by h 
(step S324). 

Through the steps described above, it follows that the lith local surface 
corresponds to the current source closest to the brain surface, the initial estimated value 
of the depth thereof is calculated as Rm(1i), and the process terminates (step S326). 

(II-3) Identification of depth of current source corresponding to each local 

surface 

Figs. 7 and 8 are flow charts representing a process of step SI 08 of Fig. 4, that 
is, a process for identifying depth of a current source corresponding to each local 
surface. 

Referring to Fig. 7, the value of a variable s is set to 1 as an initial value (step 
S402). Thereafter, the variable s is compared with a possible maximum value L of 
variable s (step S404), and when the variable s is not larger than the maximum value L, 
the process proceeds to step S406, and when the variable s exceeds the maximum value 
L, the process proceeds to step S434 (step S404). Specifically, the process from step 
S406 to step S432 is repeated from s=l to s=L. 

In step S406, when the depths of local surfaces identified so far are represented 
as {Rm(U), Rm(1s)}, the depths of these local surfaces are fixed at {RmOO, . Rm(1 s )}, 
respectively (step S406). 

First, it is assumed that 1 is not equal to any of {U,..., l s }, and that 1 belongs to a 
set of {1, L} (step S408). Then, whether all possible values are considered as the 
value of variable 1 or not is determined (step S410), and if all the possible values have 
been considered, the process proceeds to step S430. Otherwise, the process proceeds 
to step S412. 

In step S412, the depth of a local surface that is different from the 1th local 
surface and not any of {U, l s } is fixed at R M . 
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The value of variable k is set to 1 (step S414). Thereafter, the variable k is 
compared with a possible maximum value K of variable k (step S416), and when the 
variable k is not larger than the maximum value K, the process proceeds to step S418, 
and when the variable k exceeds the maximum value K, the process proceeds to step 
S426. Specifically, the process from step S41 8 to step S424 is repeated from k=l to 
k=K. 

In step S418, the depth of the 1th local surface is set to Rk. 

Thereafter, for the set of L local surfaces, current distribution is estimated (step 

S420). 

Further, posterior probability (free energy) for this arrangement of the local 
surfaces is calculated (step 422). 

Referring to Fig. 8, in step S426, after the process to k=K is finished in the 
above-described manner, the depth R M (1) of the 1th local surface attaining the highest 
posterior probability is calculated (step S426). 

Then, other variable 1 that is not equal to any of {U, l s }> belongs to the set of 
{ 1, L} and is not yet processed is selected (step S428), and the process returns to 
step S410. 

Through the above-described steps, among R M (1) values corresponding to all 
the processed variables 1, the value 1 that is closest to the brain surface is calculated and 
denoted by l s +i (step S430). 

Further, s is replaced by s+1 (step S432), and the process returns to step S404. 

When the process ends for all the possible values s, identification of the depths 
of current sources corresponding to respective local surfaces is finished (step S434). 

(II-4) Re-estimation of current source with increased resolution 

Figs. 9 and 10 are flow charts representing the process of step SI 10 shown in 
Fig. 4, that is, the process for re-estimating a position of a current source with higher 
spatial resolution. 

Referring to Figs. 9 and 10, first, numbers of local surfaces corresponding to the 
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current sources estimated using the initial resolution are re-numbered so that the 
surfaces have the numbers 1, 2, 3, L starting from the one closest to the brain surface 
(step S502). Then, corresponding to the expansion of current distribution on each 
local surface, the local surface is made smaller (step S504). 

The depth of the local surface RmO) calculated by the process up to step SI 08 of 
Fig. 4 is denoted as R M ° ld (l) (step S506). 

New resolution and search width in the depth direction are respectively 
represented as AR and (K L » AR). Further, resolution of lattice points on the local 
surface is also increased (step S508). 

The value of a variable 1 is set to 1 as an initial value (step S510). Thereafter, 
the variable 1 is compared with a possible maximum value L of variable 1 (step S512), 
and when the variable 1 is not larger than the maximum value L, the process proceeds to 
step S514, and when the variable 1 exceeds the maximum value L, the process proceeds 
to step S534 (step S512). Specifically, the process from step S514 to step S532 is 
repeated from 1=1 to 1=L. 

In step S514, the depth of a local surface 1' other than the 1th local surface is 
fixed to R M old (1'). 

Further, the process from step S518 to step S526 below is performed with k=0, 
±1,...,±K L . 

First, the value of variable k is set to 0 (step S516). Thereafter, the absolute 
value of variable k is compared with the possible maximum absolute value of K L of 
variable k (step S518), and when the absolute value of variable k is not larger than the 
maximum value K L , the process proceeds to step S520, and when the absolute value of 
variable k exceeds the maximum value K L , the process proceeds to step S528 (step 
S518). 

In step S520, the depth of the 1th local surface is set to (R M old (l)+k»AR). 
Thereafter, for the set of L local surfaces, current distribution is estimated using 
the present resolution (step S522). 
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Further, posterior probability (free energy) for this arrangement of the local 
surfaces is calculated (step S524). Then, the value k is set to the next one of {±1, 
±K L }, and the process returns to step S518. 

After the above-described process is performed until k=±K L , the value k that 
provides the highest posterior probability is calculated in step S528, which value is 
denoted by k M . 

Then, R M ° ld (l) is replaced by R M old (l)+k M *AR (step S530). The value of 
variable 1 is incremented by only 1 (step S532), and the process returns to step S512. 

When the above-described process has been performed until the value of 
variable 1 attains to the maximum value L and the resolution has reached the final 
resolution (step S534), the process is terminated (step S538). 

When the resolution is not yet the final resolution (step S534), the resolution of 
the lattice points on the local surface and the resolution in the depth direction are 
increased (step S536), and the process returns to step S510. 

By the "method of estimating brain current source" as described above, it 
becomes possible to estimate the position of a brain current source, including the 
position in the depth direction, from observation data of MEG (or EEG). Further, 
such estimation in the depth direction is applicable even when there are a plurality of 
current sources. Still further, the method is applicable when the current sources are 
localized as in the case of current dipoles or when the current source has an expansion. 
In addition, by the method, it is also possible to estimate how the current source 
expands. 

By additionally performing the process described with reference to Figs. 9 and 
10 after the estimation with the initial resolution, it becomes possible to successively 
increase resolution for estimating a position. This also means that it is possible to 
review with the scope of search restricted based on physiological findings or the like. 

[Result of simulation] 

In the following, simulation results will be described, in which current source 
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positions of an assumed model were estimated in accordance with the method of 
estimating brain current source described above. 

Fig. 1 1 is a top view of a magnetic field distribution in the radial direction 
observed on a surface of a hemisphere, assuming that the human brain is a hemisphere 
having the radius of 10.0 (arbitrary unit). 

Fig. 1 1 shows a magnetic field distribution on a surface of the hemisphere in 
which a single current source is positioned at a radius of 7.0 from the center, as will be 
described later. In the simulation that will be discussed below, noise having the S/N 
ratio of 0, 1 is added to the observation data of the magnetic field. 

Fig. 12 represents results of initial estimation of current sources using initial 
resolution. 

Fig. 12 shows the result of calculation of current distribution, in which, as 
described with reference to the process of steps SI 02 to SI 04 of Fig. 4, in order to 
perform initial estimation of the current source using the initial resolution, the depth 
(radius) is changed stepwise and on the hemisphere of each depth, the current 
distribution that attains the maximum free energy is calculated in accordance with 
variational Bayes method. 

Fig. 12(a) shows the result where the radius is R=5.0, 12(b) shows the result 
where the radius is R=6.0 and 12(c) shows the result where the radius is R=7.0. 

Fig. 13 represents results of initial estimation of current sources using initial 
resolution, where the radius is larger (that is, closer to the surface of the brain). 

Fig. 13(d) shows the result where the radius is R=7.5, 13(e) shows the result 
where the radius is R=8.0 and 13(f) shows the result where the radius is R=9.0. 

Fig. 14 represents magnitude of free energy on each virtual hemispherical 
surface, obtained when current distribution that attains maximum free energy is 
calculated for virtual hemispherical surfaces at various depths assumed in the brain. 

Referring to Fig. 14, it can be understood that the maximum point of free energy 
exists between the radii of 7 and 8, when calculation is made assuming that the 
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hemispherical surface as a whole is a virtual curved surface. 

As can be seen from Fig. 13(d), one point of relative maximum exists in the 
current distribution on the virtual hemispherical surface having the radius of R=7.5. 

As the current source is initially estimated using the initial resolution in this 
manner, on the virtual curved surface having the maximum free energy, a local surface 
that encompasses the relative maximum point of current distribution is calculated, and 
the process for identifying the current source is performed on the local surface. 

Figs. 15 and 16 represent processes for identifying the current source executed 
further on a local surface including a point of relative maximum of the current 
distribution, with the resolution of lattice points on the local surface and the resolution 
in the depth direction increased. 

Figs. 15(a), 15(b) and 15(c) represent current distribution on local surfaces 
where the radius is R=5.0, R=6.0 and R=7.0, respectively. Similarly, Figs. 16(d), 
16(e) and 16(f) represent current distribution where the radius is R=7.5, R=8.0 and 
R=9.0, respectively. 

Fig. 17 represents the resultant free energy calculated with the depth of local 
surface varied. 

Referring to Fig. 17, the free energy as a result of calculation using the local 
surface has the relative maximum value near the radius of R=7, and from the result, it is 
possible to identify that one current source exists at the position of R=7. 

The current distribution of the current source at this time is as shown in Fig. 
15(c), and as already described, it can be understood that the current distribution also 
has the narrowest expansion. 

Fig. 18 is a top view of a magnetic field distribution observed on a surface of 
the brain assumed as a hemisphere, when there are two current sources in the brain. 

In the example shown in Fig. 18, it is assumed that the current sources are 
positioned at the radii of R=6.0 and R=8.0. 

Figs. 19 and 20 represent results of initial estimation of current sources using 
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initial resolution. 

Specifically, Figs. 19 and 20 illustrate the process (step SI 04) for performing 
initial estimation of current sources using the initial resolution shown in Fig. 4, when 
there are two current sources in the brain. 

Figs. 19(a), 19(b) and 19(c) represent current distribution on virtual 
hemispherical surfaces where the radius is R=5.0, R=6.0 and R=7.0, respectively. 

Figs. 20(d), 20(e) and 20(f) represent current distribution on virtual 
hemispherical surfaces where the radius is R=7.5, R=8.0 and R=9.0, respectively. 

Fig. 21 represents radius dependency of free energy, when current distribution 
that attains maximum free energy is calculated for each virtual hemispherical surface of 
the whole hemispherical surfaces. 

As shown in Fig. 21, when the free energy is calculated over the entire 
hemispherical surface, it is understood that the maximum value of free energy exists 
between the radii of R=7 and R=8, that is, near the radius of R=7.5. 

As can be seen from Fig. 20(d), on the virtual hemispherical surface having the 
radius of R=7.5, there are two relative maximum points in the current distribution and 
two local surfaces are determined. 

Next, the result of calculation corresponding to the process of step SI 06 shown 
in Fig. 4 will be described. 

Figs. 22 and 23 illustrate the process for calculating the depth of an 1th local 
surface R M (1) attaining maximum posterior probability, to identify a current source 
closest to the brain surface. 

Therefore, Figs. 22 and 23 represent current distribution when a local surface 
corresponding to one current source is moved while the other local surface is fixed on a 
surface that attains the maximum free energy calculated over the entire hemispherical 
surface. 

Figs. 22(a), 22(b) and 22(c) represent current distribution on local surfaces 
where the radius is R=5.0, R=6.0 and R=7.0, respectively. 
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Figs. 23(d), 23(e) and 23(f) represent current distribution on local surfaces 
where the radius is R=7.5, R=8.0 and R=9.0, respectively. 

Fig. 24 represents local surface position (radius) dependency of free energy 
when a virtual local surface is moved as shown in Figs. 22 and 23. 

As shown in Fig. 24, when a local surface corresponding to one current source 
is moved to the position of radius R=8 while the other local surface is fixed on a surface 
that attains the maximum free energy calculated over the entire hemispherical surface, 
the free energy is maximized. 

Therefore, from the result shown in Fig. 24, it can be seen that the position of 
the current source closest to the brain surface is at the radius of R=8.0. 

Thereafter, the depth of the other current source is identified. 

Here, the depth of the local surface corresponding to one current source is fixed 
at the radius of R=8.0, and the depth of the local surface corresponding to the other 
current source is varied. 

Figs. 25 and 26 represent current distribution when one local surface is fixed on 
a specific surface and a local surface corresponding to the other current source is moved. 

Figs. 25(a), 25(b) and 25(c) represent current distribution on local surfaces 
where the radius is R=5.0, R=5.5 and R=6.0, respectively. 

Figs. 26(d), 26(e) and 26(f) represent current distribution on local surfaces 
where the radius is R=6.5, R=7.0 and R=7.5, respectively. 

Fig. 27 represents local surface position (radius) dependency of free energy 
when a virtual local surface is moved as shown in Figs. 25 and 26. 

As can be seen from Fig. 27, here, the free energy attains relative maximum at 
the position of radius R=6. 

It can be understood that by such a process, even when there are two current 
sources in the brain, it is possible to identify the depth of each of the current sources. 

It should be understood that the embodiments disclosed herein are illustrative 
and non-restrictive in every respect. The scope of the present invention is defined by 
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the terms of the claims, rather than the description above, and is intended to include any 
modifications within the scope and meaning equivalent to the terms of the claims. 
[Effects of the Invention] 

As described above, according to the brain current source estimating method of 
the present invention, it is possible to estimate the position of a brain current source 
including the depth direction, from observation data of MEG or EEG. Such 
estimation in the depth direction is applicable even when there are a plurality of current 
sources. Further, the method is applicable where the current sources exist locally as in 
the case of current dipoles and applicable to a current source that has a certain 
expansion. Further, it is possible to estimate how the current source expands. In 
addition, after the estimation with the initial resolution, it is possible to successively 
increase resolution for estimating a position. 
[Brief Description of the Drawings] 

Fig. 1 is a schematic illustration of an exemplary configuration of an MEG 

system. 

Fig. 2 is a schematic illustration representing an electromagnetic field generated 
by a current source, observed on an appropriate curved surface. 

Fig. 3 is a schematic illustration representing a relation between a current source 
in the brain and a plurality of virtual curved surfaces. 

Fig. 4 is a flow chart representing an overall procedure of the brain current 
source estimating method in accordance with the present invention. 

Fig. 5 is a flow chart representing a process of initial estimation of the current 
source using an initial resolution. 

Fig. 6 is a flow chart representing a process for identifying a current source 
closest to the brain surface. 

Fig. 7 is a first flow chart representing a process for identifying depth of a 
current source corresponding to each local surface. 

Fig. 8 is a second flow chart representing a process for identifying depth of a 
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current source corresponding to each local surface. 

Fig. 9 is a first flow chart representing a process for re-estimating a position of a 
current source with higher spatial resolution. 

Fig. 10 is a second flow chart representing a process for re-estimating a position 
of a current source with higher spatial resolution. 

Fig. 1 1 is a top view of a magnetic field distribution observed on a surface of a 
hemisphere, assuming that the human brain is a hemisphere. 

Fig. 12 represents results of initial estimation of current sources using initial 
resolution. 

Fig. 13 represents results of initial estimation of current sources using initial 
resolution, where the radius is greater (closer to the surface of the brain). 

Fig. 14 represents magnitude of free energy on each virtual hemispherical 
surface, obtained when current distribution that attains highest free energy is calculated. 

Fig. 1 5 is a first diagram representing processes for identifying the current 
source executed further on a local surface including a point of relative maximum. 

Fig. 16 is a second diagram representing processes for identifying the current 
source executed further on a local surface including a point of relative maximum. 

Fig. 17 represents the resultant free energy calculated with the depth of local 
surface varied. 

Fig. 18 is a top view of a magnetic field distribution observed on a surface of 
the brain assumed as a hemisphere, when there are two current sources in the brain. 

Fig. 19 is a first diagram representing results of initial estimation of current 
sources using initial resolution. 

Fig. 20 is a second diagram representing results of initial estimation of current 
sources using initial resolution. 

Fig. 21 represents radius dependency of free energy, when current distribution 
that attains highest free energy is calculated for each virtual hemispherical surface. 

Fig. 22 is a first diagram for illustrating the process for calculating the depth of 
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an 1th local surface attaining maximum posterior probability, to identify a current 
source closest to the brain surface. 

Fig. 23 is a second diagram for illustrating the process for calculating the depth 
of an 1th local surface attaining maximum posterior probability, to identify a current 
source closest to the brain surface. 

Fig. 24 represents local surface position (radius) dependency of free energy 
when a virtual local surface is moved. 

Fig. 25 is a first diagram representing current distribution when one local 
surface is fixed on a specific surface and a local surface corresponding to the other 
current source is moved. 

Fig. 26 is a second diagram representing current distribution when one local 
surface is fixed on a specific surface and a local surface corresponding to the other 
current source is moved. 

Fig. 27 represents local surface position (radius) dependency of free energy 
when a virtual local surface is moved. 

[Description of the Reference Characters] 

10 subject; 12 MEG; and 20 computer system. 



-43 - 



JP 2001-400519 




FIG.4 



ESTIMATION OF CURRENT SOURCE 





START 


S100 




\ 


f 


DETERMINE PARAMETER 
DESIGNATING SHAPE OF HIERARCHICAL 
PRIOR DISTRIBUTION 


\ 


r 


INITIAL ESTIMATION OF CURRENT 
SOURCE USING INITIAL RESOLUTION 
(SEARCH FOR LOCALIZED CURRENT 
SOURCE CANDIDATES) 


\ 


f 


IDENTIFY CURRENT SOURCE 
CLOSEST TO BRAIN SURFACE 


\ 




IDENTIFY DEPTH OF OTHER CURRENT SOURCE 


\ 


f 


RE-ESTIMATION 

WITH HIGHER SPATIAL RESOLUTION 




\ 


i 






END 


S112 



S102 



S104 



S106 



S108 



S110 



FIG.5 



INITIAL ESTIMATION OF 
CURRENT SOURCE 



START 

(INITIAL ESTIMATION OF 
CURRENT SOURCE) 



ESTIMATE CURRENT DISTRIBUTION FOR VIRTUAL CURVED 
SURFACES OF RESPECTIVE DEPTHS OF SAMPLE POINTS 
{Rk} IN DEPTH DIRECTION WITH INITIAL SPATIAL RESOLUTION 

OBTAIN POSTERIOR PROBABILITY (FREE ENERGY) 
FOR EACH DEPTH Rk 


S202 


\ 


1 






CALCULATE POINT OF RELATIVE MAXIMUM OF CURRENT 
INTENSITY FOR CURRENT DISTRIBUTION AT DEPTH Rm 
THAT ATTAIN HIGHEST POSTERIOR PROBABILITY 


S204 




l 






DETERMINE LOCAL SURFACE INCLUDING EACH POINT 

OF RELATIVE MAXIMUM (NUMBER OF LOCAL SURFACES: L) 


S206 






1 








END 


S208 





FIG.6 



START 

(ESTIMATION OF CURRENT SOURCE 
CLOSEST TO BRAIN SURFACE) 

i 



ESTIMATION OF CURRENT 
SOURCE CLOSEST TO 
BRAIN SURFACE 



S300 




S302 



S304 



FIX DEPTH OF LOCAL SURFACES 
OTHER THAN ITH LOCAL SURFACE TO Rm 




S306 



S308 



S310 



SET DEPTH OF ITH LOCAL SURFACE TO Rk 



I 



PERFORM CURRENT DISTRIBUTION ESTIMATION 
ON A SET OF L LOCAL SURFACES 



T 



CALCULATE POSTERIOR PROBABILITY FOR 
THIS ARRANGEMENT OF LOCAL SURFACES 



± 



S312 



S314 



S316 



k=k+1 



S318 



CALCULATE DEPTH Rm(I) OF ITH LOCAL SURFACE 
THAT ATTAINS HIGHEST POSTERIOR PROBABILITY 



S320 



1=1+1 



S322 



CALCULATE l(=H) CLOSEST TO BRAIN SURFACE AMONG Rm(I) 



I 



S324 



END 



JS326 



FIG.7 



START 

(IDENTIFICATION OF DEPTH OF 
CURRENT SOURCE) 

* : 



IDENTIFICATION OF DEPTH 
OF CURRENT SOURCE 



S400 




S402 



S404 



FIX ALREADY IDENTIFIED DEPTHS {Rm(Ii),...,Rm(Is)} 
OF LOCAL SURFACES TO RESPECTIVE VALUES 



S406 



Is} AND I e{1 L} 




S408 



S410 



FIX DEPTH OF LOCAL SURFACE OTHER THAN 
ITH CURVED SURFACE AND NOT {h Is} TO Rm 



S412 




S414 



S416 



SET DEPTH OF ITH LOCAL SURFACE TO Rk 



I 



S418 



PERFORM CURRENT DISTRIBUTION ESTIMATION 
FOR THE SET OF L LOCAL SURFACES 



I 



S420 



CALCULATE POSTERIOR PROBABILITY FOR 
THIS ARRANGEMENT OF LOCAL SURFACES 



S422 



k=k+1 
I — 



S424 



FROM S428 
FROM S432 



V 

TO S426 



TO S430 

TO S434 



FIG.8 



TO S410 



TO S404 
A 



FROM S410 

FROM S404 



CALCULATE DEPTH Rm(I) OF ITH LOCAL SURFACE 
THAT ATTAINS HIGHEST POSTERIOR PROBABILITY 



S426 



SELECT OTHER I OF l^{h,...,| s } AND I L} 

I 



S428 



CALCULATE l(=ls+l) CLOSEST 
TO BRAIN SURFACE AMONG Rm(I) 



S430 



s=s+1 



S432 



END 



S434 



FIG.9 



RE-ESTIMATION OF 
CURRENT SOURCE WITH 
HIGHER RESOLUTION 



START 

(RE-ESTIMATION OF CURRENT 
SOURCE WITH HIGHER RESOLUTION) 



S500 



RENUMBER LOCAL SURFACES CORRESPONDING TO 
CURRENT SOURCES ESTIMATED WITH INITIAL RESOLUTION, 
STARTING FROM ONE CLOSEST TO BRAIN SURFACE 



S502 



REDUCE LOCAL SURFACE IN CORRESPONDENCE TO EXPIATION 
OF CURRENT DISTRIBUTION OF EACH LOCAL SURFACE 



S504 



SET DEPTH Rm(I) ALREADY IDENTIFIED 
WITH INITIAL RESOLUTION TO DEPTH RM old (l) 



J. 



S506 



INCREASE RESOLUTION OF LOCAL SURFACE, SET NEW DEPTH 
DIRECTION RESOLUTION AND SEARCH WIDTH TO AR AND (Kl- AR) 



S508 



c 



1=1 




]S510 



S512 



FIX DEPTH OF LOCAL SURFACE 

I' OTHER THAN ITH LOCAL SURFACE TO RM old (l') 



FROM 
S536 




S514 



S518 



SET DEPTH OF ITH LOCAL SURF ACE TO (RM old (l)+k - A R) 

i 



S520 



PERFORM CURRENT DISTRIBUTION ESTIMATION ON A 
SET OF L LOCAL SURFACES WITH PRESENT RESOLUTION 

j S522 



CALCULATE POSTERIOR PROBABILITY FOR 
THIS ARRANGEMENT OF LOCAL SURFACES 

j 



S524 



SET k TO NEXT VALUE OF {±1,..., ±Kl} 



S526 



FROM 
S532 



TO S528 



TO 
S534 



FIG.10 



TO S510 TO S512 



FROM 



S512 



FROM S518 



A 



CALCULATE k THAT ATTAINS HIGHEST 
POSTERIOR PROBABILITY AND SET THIS TO kM 


\ 


I 


S528 


RM old (l)=RM old (l)+kM- AR 




\ 


1 


S530 




1=1+1 


S532 



RESOLUTION 
IS FINAL RESOLUTION 

9 

~S534 

N 



INCREASE LOCAL SURFACE RESOLUTION 
AND DEPTH DIRECTION RESOLUTION A R 



S536 



END 



S538 



FIG.11 




FIG.12(a) 



R=5.0 


( 




FIG.12(b) 


R=6.0 




FIG.12(c) 


R=7.0 











FlG.13(d) 



R=7.5 




FIG. 13(e) 


R=8.0 






FIG. 13(f) 


R=9.0 







FIG. 15(a) 



R=5.0 



FIG. 15(b) 



R=6.0 



FIG. 15(c) 



R=7.0 



FIG. 16(d) 




FIG. 16(f) 



R=9.0 




FIG. 18 




FIG. 19(a) 




FIG.20(d) 




FIG.21 




FIG.22(a) 



R=5.0 




FIG.22(b) 



R=6.0 




FIG.22(c) 



FIG.23(d) 



R=7.5 




FIG.23(e) 




FIG.24 




FIG.25(a) 



R=5.0 



FIG.25(b) 



R=5.5 



FIG.25(c) 



R=6.0 



FIG.26(d) 



R=6.5 




FIG.26(e) 




JP2001-400519 



[Document Name] Abstract 
[Abstract] 

[Subject] A brain current source estimating method that enables estimation of a 
position or distribution, depth direction inclusive, of a brain current source from 
observed data is provided. 

[Solving Means] Based on the fact that by causing an appropriate current flow on a 
virtual curved surface between a current source and an observation surface, an 
electromagnetic field generated by the current source can be recovered, and that as the 
curved surface comes closer to the true current source, expansion of current distribution 
on the curved surface becomes smaller, Bayesian estimation of the current source that 
recovers the observed data is performed. In this estimation, the fact that the model 
posterior probability attains the maximum when the curved surface includes the current 
source is utilized, that is, the model posterior probability is considered, to estimate the 
position of the current source including the depth direction. 
[Selected Drawing] Fig. 3 
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